To record the the co-occurence network of pest injuries by country
library(ggplot2)
require(dplyr)
## Loading required package: dplyr
##
## Attaching package: 'dplyr'
##
## The following objects are masked from 'package:stats':
##
## filter, lag
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
require(reshape)
## Loading required package: reshape
##
## Attaching package: 'reshape'
##
## The following object is masked from 'package:dplyr':
##
## rename
require(reshape2)
## Loading required package: reshape2
##
## Attaching package: 'reshape2'
##
## The following objects are masked from 'package:reshape':
##
## colsplit, melt, recast
require(qgraph)
## Loading required package: qgraph
library(gridExtra)
library(cowplot)
##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:ggplot2':
##
## ggsave
library(lubridate)
##
## Attaching package: 'lubridate'
##
## The following object is masked from 'package:reshape':
##
## stamp
library(XLConnect)
## Loading required package: XLConnectJars
## XLConnect 0.2-11 by Mirai Solutions GmbH [aut],
## Martin Studer [cre],
## The Apache Software Foundation [ctb, cph] (Apache POI, Apache Commons
## Codec),
## Stephen Colebourne [ctb, cph] (Joda-Time Java library)
## http://www.mirai-solutions.com ,
## http://miraisolutions.wordpress.com
You can also embed plots, for example:
Filepath <- "E:/Google Drive/1.SKEP1/SKEP1survey.xls" # please check your file in shared google drive
Filepath <- "~/Google Drive/1.SKEP1/SKEP1survey.xls"
data <- readWorksheetFromFile(Filepath, sheet = 1)
#### clean define the missing value ####
data[data == "-"] <- NA # replace '-' with NA
data[data == ""] <- NA # replace 'missing data' with NA
#### end cleaning of data ####
#### to lower variable names ####
names(data) <- tolower(names(data))
#### end setting the varibales ####
data <- transform(data,
phase = as.factor(phase),
fno = as.character(fno),
identifier = as.character(identifier),
country = as.factor(country),
year = as.factor(year),
season = as.factor(season),
lat = as.numeric(lat),
long = as.numeric(long),
village = as.character(village),
fa = as.numeric(fa),
fn = as.character(fn),
lfm = as.character(lfm),
pc = as.factor(pc),
fp = as.character(fp),
cem = as.factor(cem),
ast = as.factor(ast),
nplsqm = as.numeric(nplsqm),
ced = dmy(ced),# Date data try to use as.Data(., format = '%d-%b-%y') it is not working
cedjul = as.numeric(cedjul),
hd = dmy(hd),
hdjul = as.numeric(hdjul),
ccd = as.numeric(ccd),
cvr = as.character(cvr),
vartype = as.factor(vartype),
varcoded = as.factor(varcoded),
fym = as.character(fym),
fymcoded = as.factor(fymcoded),
n = as.numeric(n),
p = as.numeric(p) ,
k = as.numeric(k),
mf = as.numeric(mf),
wcp = as.factor(wcp),
mu = as.character(mu) ,
iu = as.numeric(iu),
hu = as.numeric(hu),
fu = as.numeric(fu),
cs = as.factor(cs),
ldg = as.numeric(ldg),
yield = as.numeric(yield) ,
dscum = as.factor(dscum),
wecum = as.factor(wecum),
ntmax = as.numeric(ntmax),
npmax = as.numeric(npmax),
nltmax = as.numeric(nltmax),
nlhmax = as.numeric(nltmax),
waa = as.numeric(waa),
wba = as.numeric(wba) ,
dhx = as.numeric(dhx),
whx = as.numeric(whx),
ssx = as.numeric(ssx),
wma = as.numeric(wma),
lfa = as.numeric(lfa),
lma = as.numeric(lma),
rha = as.numeric(rha) ,
thrx = as.numeric(thrx),
pmx = as.numeric(pmx),
defa = as.numeric(defa) ,
bphx = as.numeric(bphx),
wbpx = as.numeric(wbpx),
awx = as.numeric(awx),
rbx =as.numeric(rbx),
rbbx = as.numeric(rbbx),
glhx = as.numeric(glhx),
stbx=as.numeric(stbx),
rbpx = as.numeric(rbpx),
hbx= as.numeric(hbx),
bbx = as.numeric(bbx),
blba = as.numeric(blba),
lba = as.numeric(lba),
bsa = as.numeric(bsa),
blsa = as.numeric(blsa),
nbsa = as.numeric(nbsa),
rsa = as.numeric(rsa),
lsa = as.numeric(lsa),
shbx = as.numeric(shbx) ,
shrx = as.numeric(shrx),
srx= as.numeric(srx),
fsmx = as.numeric(fsmx),
nbx = as.numeric(nbx),
dpx = as.numeric(dpx),
rtdx = as.numeric(rtdx),
rsdx = as.numeric(rsdx),
gsdx =as.numeric(gsdx),
rtx = as.numeric(rtx)
)
## Warning: All formats failed to parse. No formats found.
## Warning: NAs introduced by coercion
#### Delete the unnessary variables variables without data (NA) ####
data$phase <- NULL # there is only one type yype of phase in the survey
data$identifier <- NULL # this variable is not included in the analysis
data$village <- NULL
data$fa <- NULL # field area is not include in the analysis
data$fn <- NULL # farmer name can not be included in this survey analysis
data$fp <- NULL # I do not know what is fp
data$lfm <- NULL # there is only one type of land form in this survey
data$ced <- NULL # Date data can not be included in the network analysis
data$cedjul <- NULL
data$hd <- NULL # Date data can not be included in the network analysis
data$hdjul <- NULL
data$cvr <- NULL
data$varcoded <- NULL # I will recode them
data$fym.coded <- NULL
data$mu <- NULL # no record
data$nplsqm <- NULL
#### Delete the unnessary variables variables without data (NA) ####
OutputProfile <- data %>%
select (
country,
season,
dhx,
whx,
ssx,
wma,
lfa,
lma,
rha,
thrx,
pmx,
defa,
bphx,
wbpx,
awx,
rbx,
rbbx,
glhx,
stbx,
hbx,
bbx,
blba,
lba,
bsa,
blsa,
nbsa,
rsa,
lsa,
shbx,
shrx,
srx,
fsmx,
nbx,
dpx,
rtdx,
rsdx,
gsdx,
rtx
)
data <- OutputProfile
#### The Philippines n = 40 ####
php <- data %>%
filter( country == "PHL") %>%
select(-c(country, season))
php <- php[,apply(php, 2, var, na.rm = TRUE) != 0] # exclude the column with variation = 0
php <- php[complete.cases(php),] # exclude row which cantain NA
#### end Philippines subset ####
#### India N = 105 ####
ind <- data %>%
filter( country == "IND") %>%
select(-c(country, season))
ind <- ind[,apply(ind, 2, var, na.rm = TRUE) != 0] # exclude the column with variation = 0
ind <- ind[complete.cases(ind),] # exclude row which cantain NA
#### end India subset ####
##### Indonesia N= 100 ####
idn <- data %>%
filter( country == "IDN") %>%
select(-c(country, season))
idn <- idn[,apply(idn, 2, var, na.rm = TRUE) != 0] # exclude the column with variation = 0
idn <- idn[complete.cases(idn),] # exclude row which cantain NA
#### end Indonesia subset ####
#### Thailand n = 105 ####
tha <- data %>%
filter( country == "THA") %>%
select(-c(country, season))
tha <- tha[,apply(tha, 2, var, na.rm = TRUE) != 0] # exclude the column with variation = 0
tha <- tha[complete.cases(tha),] # exclude row which cantain NA
#### end Thailand subset ####
#### Vietnam n = 105 ####
vnm <- data %>%
filter( country == "VNM") %>%
select(-c(country, season))
vnm <- vnm[,apply(vnm, 2, var, na.rm = TRUE) != 0] # exclude the column with variation = 0
vnm <- vnm[complete.cases(vnm),]
abbre <- c("DH", "WH", "SS", "WM", "LF", "DEF", "BPH", "WPH", "AM", "RB", "GLH", "BLB", "LB", "BS", "BLS", "NBS", "RS", "LS", "SHB", "SHR", "SR", "FSM", "NB", "DP", "RTD", "RT")
php$ldg <- NULL
names(php) <-abbre
InsectInjuiry <- c(1:11, 26)
Disease <- 12:25
injuries <- list(InsectInjuiry, Disease)
#==========#
cormat <- cor(php, method = "spearman", use = "pairwise") # spearman correlation
php_qnet <- qgraph(cormat,
layout = "spring",
threshold = 0.20,
maximum = 1,
groups = injuries,
color = c("skyblue", "wheat"),
vsize = 5,
line = 3,
diag = FALSE,
posCol = "forestgreen",
negCol = "firebrick3",
borders = FALSE,
# legend = TRUE,
vTrans = 200,
#nodeNames = Names,
legend.cex = 0.3,
title = "Philippines"
#filetype = "png",
#filenames = "figs/ind_network.png"
)
qgraph(cormat,
layout = "spring",
threshold = 0.20,
maximum = 1,
groups = injuries,
color = c("skyblue", "wheat"),
vsize = c(1.5,10),
line = 3,
diag = FALSE,
posCol = "forestgreen",
negCol = "firebrick3",
borders = FALSE,
# legend = TRUE,
vTrans = 200,
#nodeNames = Names,
legend.cex = 0.3,
title = "Philippines"
#filetype = "png",
#filenames = "figs/ind_network.png"
)
# compute the topological properties
cen <- centrality_auto(php_qnet)
topo <- cen$node.centrality
topo$node <- row.names(topo)
row.names(topo) <- NULL
topo.melt <- melt(topo)
## Using node as id variables
p1 <- topo %>% ggplot(aes(x= Strength, y = reorder(node, Strength))) +
geom_point(size = 3, color ="red") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Node degree") +
ylab("Variables")
p2 <- topo %>% ggplot(aes(x= Betweenness, y = reorder(node, Betweenness))) +
geom_point(size = 3, color ="blue") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Betweenness") +
ylab("Variables")
# Compute clustering coefficients
clust <- clustcoef_auto(php_qnet)
clusZ.sp <- NULL
clusZ.sp$node <- rownames(clust)
clusZ.sp$clustZhang <- clust$clustZhang
clusZ.sp <- as.data.frame(clusZ.sp)
p3 <- clusZ.sp %>%
ggplot(aes(x= clustZhang, y = reorder(node, clustZhang))) +
geom_point(size = 3, color ="yellow") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Clustering coefficient") +
ylab("Variables")
plot_grid(p1, p2, p3, labels=c("A", "B", "C"), ncol = 3, nrow = 1 )
abbre <- c("DH", "WH", "SS", "WM", "LF", "DEF", "BPH", "WPH", "AM", "RB", "RBB", "GLH", "STB", "BLB", "LB", "BS", "BLS", "NBS", "RS", "LS", "SHB", "SHR", "SR", "FSM", "NB", "DP", "RTD", "RSD", "RT")
names(idn) <- abbre
InsectInjuiry <- c(1:13, 29)
Disease <- 14:28
injuries <- list(InsectInjuiry, Disease)
#==========#
cormat <- cor(idn, method = "spearman", use = "pairwise") # spearman correlation
idn_qnet <- qgraph(cormat,
layout = "spring",
threshold = 0.20,
maximum = 1,
groups = injuries,
color = c("skyblue", "wheat"),
vsize = 5,
line = 3,
diag = FALSE,
posCol = "forestgreen",
negCol = "firebrick3",
borders = FALSE,
# legend = TRUE,
vTrans = 200,
#nodeNames = Names,
legend.cex = 0.3,
title = "Indonesia"
#filetype = "png",
#filenames = "figs/ind_network.png"
)
qgraph(cormat,
layout = "spring",
threshold = 0.20,
maximum = 1,
groups = injuries,
color = c("skyblue", "wheat"),
vsize = c(1.5, 10),
line = 3,
diag = FALSE,
posCol = "forestgreen",
negCol = "firebrick3",
borders = FALSE,
# legend = TRUE,
vTrans = 200,
#nodeNames = Names,
legend.cex = 0.3,
title = "Indonesia"
#filetype = "png",
#filenames = "figs/ind_network.png"
)
# compute the topological properties
cen <- centrality_auto(idn_qnet)
topo <- cen$node.centrality
topo$node <- row.names(topo)
row.names(topo) <- NULL
topo.melt <- melt(topo)
## Using node as id variables
p1 <- topo %>% ggplot(aes(x= Strength, y = reorder(node, Strength))) +
geom_point(size = 3, color ="red") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Node degree") +
ylab("Variables")
p2 <- topo %>% ggplot(aes(x= Betweenness, y = reorder(node, Betweenness))) +
geom_point(size = 3, color ="blue") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Betweenness") +
ylab("Variables")
# Compute clustering coefficients
clust <- clustcoef_auto(idn_qnet)
clusZ.sp <- NULL
clusZ.sp$node <- rownames(clust)
clusZ.sp$clustZhang <- clust$clustZhang
clusZ.sp <- as.data.frame(clusZ.sp)
p3 <- clusZ.sp %>%
ggplot(aes(x= clustZhang, y = reorder(node, clustZhang))) +
geom_point(size = 3, color ="yellow") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Clustering coefficient") +
ylab("Variables")
plot_grid(p1, p2, p3, labels=c("A", "B", "C"), ncol = 3, nrow = 1 )
abbre <- c("DH", "WH", "WM", "LF", "DEF", "BPH", "WPH", "RB", "RBB", "GLH", "BLB", "LB", "SHB", "FSM", "NB", "RT")
names(ind) <-abbre
InsectInjuiry <-c(1:10, 16)
Disease <- 11:15
injuries <- list(InsectInjuiry, Disease)
#==========#
cormat <- cor(ind, method = "spearman", use = "pairwise") # spearman correlation
ind_qnet <- qgraph(cormat,
layout = "spring",
threshold = 0.20,
maximum = 1,
groups = injuries,
color = c("skyblue", "wheat"),
vsize = 5,
line = 3,
diag = FALSE,
posCol = "forestgreen",
negCol = "firebrick3",
borders = FALSE,
# legend = TRUE,
vTrans = 200,
#nodeNames = Names,
legend.cex = 0.3,
title = "India"
#filetype = "png",
#filenames = "figs/ind_network.png"
)
qgraph(cormat,
layout = "spring",
threshold = 0.20,
maximum = 1,
groups = injuries,
color = c("skyblue", "wheat"),
vsize = c(1.5, 10),
line = 3,
diag = FALSE,
posCol = "forestgreen",
negCol = "firebrick3",
borders = FALSE,
# legend = TRUE,
vTrans = 200,
#nodeNames = Names,
legend.cex = 0.3
#title = "Spearman's correlation"
#filetype = "png",
#filenames = "figs/ind_network.png"
)
# compute the topological properties
cen <- centrality_auto(ind_qnet)
topo <- cen$node.centrality
topo$node <- row.names(topo)
row.names(topo) <- NULL
topo.melt <- melt(topo)
## Using node as id variables
p1 <- topo %>% ggplot(aes(x= Strength, y = reorder(node, Strength))) +
geom_point(size = 3, color ="red") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Node degree") +
ylab("Variables")
p2 <- topo %>% ggplot(aes(x= Betweenness, y = reorder(node, Betweenness))) +
geom_point(size = 3, color ="blue") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Betweenness") +
ylab("Variables")
# Compute clustering coefficients
clust <- clustcoef_auto(ind_qnet)
clusZ.sp <- NULL
clusZ.sp$node <- rownames(clust)
clusZ.sp$clustZhang <- clust$clustZhang
clusZ.sp <- as.data.frame(clusZ.sp)
p3 <- clusZ.sp %>%
ggplot(aes(x= clustZhang, y = reorder(node, clustZhang))) +
geom_point(size = 3, color ="yellow") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Clustering coefficient") +
ylab("Variables")
plot_grid(p1, p2, p3, labels=c("A", "B", "C"), ncol = 3, nrow = 1 )
abbre <- c("DH", "WH", "SS", "WM", "LF", "DEF", "BPH", "WPH", "RB", "GLH", "BLB", "LB", "BS", "BLS", "NBS", "RS", "LS", "SHB", "SHR", "SR", "FSM", "NB", "DP", "RT")
names(tha) <- abbre
InsectInjuiry <- c(1:10, 24)
Disease <- 11:23
injuries <- list(InsectInjuiry, Disease)
#==========#
cormat <- cor(tha, method = "spearman", use = "pairwise") # spearman correlation
tha_qnet <- qgraph(cormat,
layout = "spring",
threshold = 0.20,
maximum = 1,
groups = injuries,
color = c("skyblue", "wheat"),
vsize = 5,
line = 3,
diag = FALSE,
posCol = "forestgreen",
negCol = "firebrick3",
borders = FALSE,
# legend = TRUE,
vTrans = 200,
#nodeNames = Names,
legend.cex = 0.3,
title = "Thailand"
#filetype = "png",
#filenames = "figs/ind_network.png"
)
qgraph(cormat,
layout = "spring",
threshold = 0.20,
maximum = 1,
groups = injuries,
color = c("skyblue", "wheat"),
vsize = c(1.5, 10),
line = 3,
diag = FALSE,
posCol = "forestgreen",
negCol = "firebrick3",
borders = FALSE,
# legend = TRUE,
vTrans = 200,
#nodeNames = Names,
legend.cex = 0.3,
title = "Thailand"
#filetype = "png",
#filenames = "figs/ind_network.png"
)
# compute the topological properties
cen <- centrality_auto(tha_qnet)
topo <- cen$node.centrality
topo$node <- row.names(topo)
row.names(topo) <- NULL
topo.melt <- melt(topo)
## Using node as id variables
p1 <- topo %>% ggplot(aes(x= Strength, y = reorder(node, Strength))) +
geom_point(size = 3, color ="red") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Node degree") +
ylab("Variables")
p2 <- topo %>% ggplot(aes(x= Betweenness, y = reorder(node, Betweenness))) +
geom_point(size = 3, color ="blue") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Betweenness") +
ylab("Variables")
# Compute clustering coefficients
clust <- clustcoef_auto(tha_qnet)
clusZ.sp <- NULL
clusZ.sp$node <- rownames(clust)
clusZ.sp$clustZhang <- clust$clustZhang
clusZ.sp <- as.data.frame(clusZ.sp)
p3 <- clusZ.sp %>%
ggplot(aes(x= clustZhang, y = reorder(node, clustZhang))) +
geom_point(size = 3, color ="yellow") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Clustering coefficient") +
ylab("Variables")
plot_grid(p1, p2, p3, labels=c("A", "B", "C"), ncol = 3, nrow = 1 )
abbre <- c("DH", "WH", "SS", "WM", "LF", "DEF", "BPH", "WPH", "AM", "RB", "GLH", "BLB", "LB", "BS", "BLS", "NBS", "RS", "LS", "SHB", "SHR", "SR", "FSM", "NB", "DP", "RTD", "RDS", "RT")
names(vnm) <-abbre
InsectInjuiry <- c(1:11, 27)
Disease <- 12:26
injuries <- list(InsectInjuiry, Disease)
#==========#
cormat <- cor(vnm, method = "spearman", use = "pairwise") # spearman correlation
vnm_qnet <- qgraph(cormat,
layout = "spring",
threshold = 0.20,
maximum = 1,
groups = injuries,
color = c("skyblue", "wheat"),
vsize = 5,
line = 3,
diag = FALSE,
posCol = "forestgreen",
negCol = "firebrick3",
borders = FALSE,
# legend = TRUE,
vTrans = 200,
#nodeNames = Names,
legend.cex = 0.3,
title = "Vietnam"
#filetype = "png",
#filenames = "figs/ind_network.png"
)
qgraph(cormat,
layout = "spring",
threshold = 0.20,
maximum = 1,
groups = injuries,
color = c("skyblue", "wheat"),
vsize = c(1.5, 10),
line = 3,
diag = FALSE,
posCol = "forestgreen",
negCol = "firebrick3",
borders = FALSE,
# legend = TRUE,
vTrans = 200,
#nodeNames = Names,
legend.cex = 0.3,
title = "Vietnam"
#filetype = "png",
#filenames = "figs/ind_network.png"
)
# compute the topological properties
cen <- centrality_auto(vnm_qnet)
topo <- cen$node.centrality
topo$node <- row.names(topo)
row.names(topo) <- NULL
topo.melt <- melt(topo)
## Using node as id variables
p1 <- topo %>% ggplot(aes(x= Strength, y = reorder(node, Strength))) +
geom_point(size = 3, color ="red") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Node degree") +
ylab("Variables")
p2 <- topo %>% ggplot(aes(x= Betweenness, y = reorder(node, Betweenness))) +
geom_point(size = 3, color ="blue") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Betweenness") +
ylab("Variables")
# Compute clustering coefficients
clust <- clustcoef_auto(vnm_qnet)
clusZ.sp <- NULL
clusZ.sp$node <- rownames(clust)
clusZ.sp$clustZhang <- clust$clustZhang
clusZ.sp <- as.data.frame(clusZ.sp)
p3 <- clusZ.sp %>%
ggplot(aes(x= clustZhang, y = reorder(node, clustZhang))) +
geom_point(size = 3, color ="yellow") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
xlab("Clustering coefficient") +
ylab("Variables")
plot_grid(p1, p2, p3, labels=c("A", "B", "C"), ncol = 3, nrow = 1 )